4(a). Time Series Decomposition

There are 3 components of time series that we have discussed:- 1. Trend. 2. Seasonality. 3. Cycles.

When we decompose a time series, we usually combine the trend and cyclic variation. Cyclic variations are difficult to seperate out because first we are going to need a lot of data to find out the actual variation that is coming from cycles. Secondly, the frequency is not fixed for the cyclic variation so even if we identify the cyclicity we cannot use for future prediction as the future cyclic might be very different from the historic one.

Time Series Component :-

An Additive Model is given by :- \[y_t = S_t\ +T_t\ + R_t \] where \(S_t\) is Seasonal Component, \(T_t\) is the trend-cycle component, and \(R_t\) is the remainder component, all at period t. Additve Decomposition is most appropriate if the magnitude of the seasonal fluctuations, or the variation around the trend-cycle, does not vary with the level of time series.

A Multiplicative Model is given by :- \[y_t = S_t\ *T_t*\ R_t.\] Additve Decomposition is most appropriate if the magnitude of the seasonal fluctuations, or the variation around the trend-cycle, appears to be proportional with the level of time series.

An alternative to multiplicative model is to transform the data to make it stable over time and then use a additive decomposition. When a log transformation has been used, this is equivalent to using a multiplicative decomposition because \[y_t = S_t\ *T_t*\ R_t => logy_t= logS_t + logT_t + logR_t.\]

Sometimes we are not interested in seasonlity in the dataset. If we remove the seasonal component from the original data, the resulting values are the “seasonally adjusted” data.

For additive Model, seasonally adjusted data is \(y_t - S_t\).

For Multiplicative Model, seasonally adjusted data is \(y_t/S_t\).

Seasonally adjusted series contains the remainder and the trend-cycle component. Therefore, they are not smooth. If the purpose is to look for turning points in a series, and interpret any changes in direction, then it is better to use the trend-cycle component rather than the seasonally adjusted data.

Moving Averages Smoothing:-

One way of decomposing a time series is to use moving average method to estimate trend-cycle. A moving average of order m can be written as \[\hat T_t = \frac{1}{m} \sum_{j = -k}^{k}y_{t+j},\ where\ m = 2k+1\] The estimate of the time-cycle component at time t is the average of k period within time t. Observations that are nearby in time are also likely to be close in value. Therefore, the average eliminates some of the randomness in the data, leaves a smooth trend-cycle component.

Example :-Annual Electricity sales data

library(fpp2)

## Loading required package: ggplot2

## Loading required package: forecast

## Loading required package: fma

## Loading required package: expsmooth

autoplot(elecsales, series = "Data") +
  autolayer(ma(elecsales, 5), series = "5-MA") +
  xlab("year") + ylab("GWh") +
  ggtitle("Annual Electricity Sales Data") +
  scale_colour_manual(values = c("Data" = "grey50", "5-MA" ="red"),
                      breaks = c("Data", "5-MA"))

## Warning: Removed 4 rows containing missing values (geom_path).

We can see that there is no estimation for the starting 2 and ending 2 years because we do not have 2 observations on the either side.

The order of the moving average decides the smoothness of the trend-cycle estimate. As we increase the order we are making the curve smoother. We usually choose the order of the MA as an odd number, so that the moving averages are symmetric.

Moving Average of moving averages :-

We can also take moving average of the moving avergaes. Someone might wants to do that to make even-order moving average symmetric.

beer2 = window(ausbeer, start = 1992)
ma4 = ma(beer2, order = 4, centre = FALSE)
ma2X4 = ma(beer2, order = 4, centre = TRUE)
ma2X4

##         Qtr1    Qtr2    Qtr3    Qtr4
## 1992      NA      NA 450.000 450.125
## 1993 450.250 446.500 446.000 443.000
## 1994 439.625 443.625 443.125 443.625
## 1995 446.125 443.875 440.375 437.000
## 1996 433.500 429.625 430.875 433.750
## 1997 434.750 438.125 440.000 439.375
## 1998 438.875 437.375 437.125 436.250
## 1999 437.125 440.250 439.000 439.625
## 2000 440.750 436.875 438.000 439.000
## 2001 436.500 435.750 431.875 432.500
## 2002 434.750 435.250 437.875 435.250
## 2003 433.625 433.500 431.500 432.750
## 2004 432.875 427.250 420.375 419.625
## 2005 420.750 423.750 430.000 430.625
## 2006 428.125 428.875 428.625 426.875
## 2007 425.125 421.500 418.375 418.375
## 2008 421.250 425.125 426.375 426.750
## 2009 428.875 430.000 429.875 426.750
## 2010      NA      NA

Moving Averages can be used to trend-cycle from seasonal data. For eg. in a 2 * 4 MA :- \[\hat T_t = \frac{1}{8}y_{t-2} + \frac{1}{4}y_{t-1} +\frac{1}{4}y_{t} +\frac{1}{4}y_{t+1} + \frac{1}{8}y_{t+2}\] When applied to quarterly data, each quarter of the year is given equal weight as first and last terms apply to the same quarter in consecutive quarter. Consequently, the seasonal variation will be averaged out and the resulting value will have little or no seasonal variation remaining.

If the seasonal period m is even, then we use 2 * m MA or (m+1) where all observations take the weight 1/m, except for the first and last terms which take weights 1/(2m) to estimate the trend-cycle variation, and if the period is odd then we can choose m-MA.

Other choices might result in trend-cycle estimate to be effected by seasonal variation.

Electrical Equipment Manufacturing :-

autoplot(elecequip, series = "Data") +
  autolayer(ma(elecequip, 12), series = "12-MA") +
  xlab("Year") + ylab("New Order Index") +
  ggtitle("Electrical Equipment Manufacturing") +
  scale_color_manual(values = c("Data" = "grey50", "12-MA" = "red"),
                     breaks = c("Data","12-MA"))

## Warning: Removed 12 rows containing missing values (geom_path).

Notice that the red line captures the trend-cycle effect without fluctuations. This shows the effectiveness of the method.

Weighted Moving Average :-

Combinations of moving averages result in weighted moving averages. For example, the 2×4-MA discussed above is equivalent to a weighted 5-MA with weights given by \([\frac{1}{8}, \frac{1}{4} , \frac{1}{4} , \frac{1}{4},\frac{1}{8}]\). In general, a weighted m-MA can be written as \[\hat T_t = \sum_{j = -k}^{k} a_j\ y_{t+j},\] where k = (m-1)/2, and a’s are the weights.It is important that the weights all sum to one and that they are symmetric so that \(a_j=a_{−j}\).

A major advantage of the weighted moving average is that the trend-cycle estimate is smooth. Instead of observations entering and leaving the calculations at full weight, their weight slowly increase or decrease as they enter and leave the calculations.

Classical Decomposition :-

In classical decomposition, we assume that the seasonal component is constant from year to year. There are 2 kind of decomposition :- Additive and Multiplicative.

1. Additive Decomposition :- 1. If m(seasonal period) is an even number, compute the trend-cycle component \(\hat T_t\) using a 2*m-MA. If m is odd, compute the trend-cycle component \(\hat T_t\) using an m-MA.

2. Calculate the de-trended series: \(y_t - \hat T_t\).

3. To estimate the seasonal component for each season, simply average the detrended values for that season. These seasonal component values are then adjusted to ensure that they add to zero. The seasonal component is obtained by stringing together these monthly values, and then replicating the sequence for each year of data. This gives \(\hat S_t\).

4. The remainder component is calculated by subtracting the estimated seasonal and trend-cycle components: \(\hat R_t = y_t − \hat T_t − \hat S_t\).

2. Multiplicative Decomposition :- 1. If m(seasonal period) is an even number, compute the trend-cycle component \(\hat T_t\) using a 2*m-MA. If m is odd, compute the trend-cycle component \(\hat T_t\) using an m-MA.

2. Calculate the de-trended series: \(y_t / \hat T_t\).

3. To estimate the seasonal component for each season, simply average the detrended values for that season. These seasonal component values are then adjusted to ensure that they add to zero. The seasonal component is obtained by stringing together these monthly values, and then replicating the sequence for each year of data. This gives \(\hat S_t\).

4. The remainder component is calculated by subtracting the estimated seasonal and trend-cycle components: \(\hat R_t = y_t / (\hat T_t \hat S_t)\).

elecequip %>% decompose(type = "multiplicative") %>%
  autoplot() + xlab("Year") +
  ggtitle("Classical Multiplicative Decomposition of Electricity Price Index")

We can see that classical approach works quite well in this case except in the year 2009 where we can see a lot of negetive remainders.

Disadvantages of Classical Decomposition :- 1. The trend-cycle estimate tends to over-smooth rapid rises and falls in the data.

2. Classical decomposition methods assume that the seasonal component repeats from year to year. For many series, this is a reasonable assumption, but for some longer series it is not.

3. The estimate of the trend-cycle is unavailable for the first few and last few observations.

4. Occasionally, the values of the time series in a small number of periods may be particularly unusual. The classical method is not robust to these kinds of unusual values.

X11 Decomposition :-

Another method to decompose a time series is the X11 method. This method overcomes the drawback of the classical method of decomposition. In particular, trend-cycle estimates are available for all observations including the end points, and the seasonal component is allowed to vary slowly over time. X11 also has some sophisticated methods for handling trading day variation, holiday effects and the effects of known predictors.

library(seasonal)
elecequip %>% 
  seas(x11 = "") -> fit
autoplot(fit) + 
  ggtitle("X11 decomposition of electrical equipment index")

We can see that the X11 method has captured the sudden fall in 2009 much better. Also, it has allowed the seasonal components to vary over the years.

Given the output from the seas() function, seasonal() will extract the seasonal component, trendcycle() will extract the trend-cycle component, remainder() will extract the remainder component, and seasadj() will compute the seasonally adjusted time series.

autoplot(elecequip, series = "Data") +
  autolayer(trendcycle(fit), series = "Trend-Cycle") +
  autolayer(seasadj(fit), series = "Seasonally Adjusted") +
  xlab("Year") + ylab("New Order Index") +
  ggtitle("Electrical equipment manufacturing") +
  scale_colour_manual(values = c("grey", "blue", "red"),
                      breaks = c("Data","Seasonally Adjusted","Trend-Cycle"))

It can be useful to use seasonal plots and seasonal sub-series plots of the seasonal component. These help us to visualise the variation in the seasonal component over time.

fit %>% seasonal() %>% ggsubseriesplot() + ylab("Seasonal")

SEATS Decomposition :-

“SEATS” stands for “Seasonal Extraction in ARIMA Time Series”. This is just another way of decomposition. The procedure works only with quarterly and monthly data.

elecequip %>% 
  seas() %>%
  autoplot() + 
  ggtitle("SEATS decomposition of electrical equipment index")

STL Decomposition :-

“Seasonal and Trend decomposition using Loess”, while Loess is a method for estimating nonlinear relationships.

STL has several advantages over the classical, SEATS and X11 decomposition methods:

1.Unlike SEATS and X11, STL will handle any type of seasonality, not only monthly and quarterly data.

2.The seasonal component is allowed to change over time, and the rate of change can be controlled by the user.

3.The smoothness of the trend-cycle can also be controlled by the user.

4.It can be robust to outliers (i.e., the user can specify a robust decomposition), so that occasional unusual observations will not affect the estimates of the trend-cycle and seasonal components. They will, however, affect the remainder component.

STL only provides facilities for additive decompositions. It is possible to obtain a multiplicative decomposition by first taking logs of the data, then back-transforming the components.

elecequip %>% 
  stl(t.window = 13, s.window = 13, robust = TRUE) %>% # periodic means that seasonal component is identical every year.
  autoplot()

Both t.window and s.window should be odd numbers; t.window is the number of consecutive observations to be used when estimating the trend-cycle; s.window is the number of consecutive years to be used in estimating each value in the seasonal component.

Measuring the strength of trend and seasonality :-

For a strongly trended data, the seasonally adjusted data should have much more variation than the remainder component. Therefore, Var(\(R_t\))/Var(\(R_t + T_t\)) should be relatively small. But for little or no trend, this variance will approximately be same.

\[F_T = max(0,1- \frac{Var(R_t)}{Var(R_t + T_t)})\] This will give a measure of the strength of the trend between 0 and 1.

The strength of seasonality is defined similarly, but with respect to the detrended data rather than the seasonally adjusted data: \[F_S = max(0,1- \frac{Var(R_t)}{Var(R_t + S_t)})\]

Forecasting with Decomposition :-

While decomposition is primarily useful for studying time series data, and exploring historical changes over time, it can also be used in forecasting.

An additive model can be written as :- \[Y_t = \hat S_t + \hat A_t\] where \(\hat A_t = \hat T_t + \hat R_t\) is the seasonally adjusted component.

Similarly, a multiplicative model can be written as :- \[Y_t = \hat S_t \hat A_t\] where \(\hat A_t = \hat T_t\hat R_t\) is the seasonally adjusted component.

We will forecast the “trend + remainder” and seasonal component seperately. It is usually assumed that the seasonal component is unchanging, or changing extremely slowly, so it is forecast by simply taking the last year of the estimated component. In other words, a seasonal naïve method is used for the seasonal component.

To forecast the seasonally adjusted component, any non-seasonal forecasting method may be used. For e.g.

fit <- stl(elecequip, t.window = 13, s.window = "periodic", robust = TRUE)
fit %>% seasadj() %>% naive() %>%
  autoplot() + xlab("New Order Index") +
  ggtitle("Naive forecasts of seasonality adjusted data")

The above graph shows the naive forecast on the seasonally adjusted data. We can reseasonalised the series by adding seasonal naive forecast.

fit %>% forecast(method = "naive") %>% #method being used on the seasonally adjusted data
  autoplot() +
  ggtitle("Naive forecasts of seasonality adjusted data")

The prediction intervals shown in this graph are constructed in the same way as the point forecasts. That is, the upper and lower limits of the prediction intervals on the seasonally adjusted data are “reseasonalised” by adding in the forecasts of the seasonal component.

A short cut approach for the same is

fcast <- stlf(elecequip, method = "naive")
autoplot(fcast)

The stlf() function uses mstl() to carry out the decomposition, so there are default values for s.window and t.window.

Excercise :-

1. Monthly sales for a product for Plastic Manufacturers

autoplot(plastics) +
  xlab("Month") + ylab("") +
  ggtitle("Sales of plastic product")

An increasing trend can be seen in the data with seasonal pattern as well. The pattern do not seems consistent for all the months.

plastics %>% decompose(type = "multiplicative") %>%
  autoplot() + xlab("Months") +
  ggtitle("Sales of plastic product")

The results above confirms what we have observed in the time plot.

plastics %>% 
  decompose(type = "multiplicative") -> fit
autoplot(seasadj(fit))

Effect of outlier on decomposition.

plastics_out = plastics
plastics_out[20] = plastics[20] + 500
plastics_out %>% 
  decompose(type = "multiplicative") -> fit
autoplot(seasadj(fit))

As we can see that the classical decomposition assumes the constant seasonal effect, the seasonality will not be effected by the outlier but the outlier can be seen in the seasonal adjusted data.

fit %>% autoplot()

As we can see that the outlier has shifted the trend a bit. Also the pattern of the seasonality has also changed a bit.

plastics_out = plastics
plastics_out[56] = plastics[56] + 500
plastics_out %>% 
  decompose(type = "multiplicative") -> fit
autoplot(seasadj(fit))

Having an outlier in the end also won’t make any difference because the seasonal variation is constant.

fit %>% autoplot()

Adding the outlier has not effected trend much as we are not giving a predition for last few values anyway. The outlier seems to have little effect on seasanol pattern as well.

2. Retail Time Series Data :-

retaildata = readxl::read_excel("/Users/atyagi/Desktop/Time Series Forecasting/Time-Series-Forecasting/Time_series_data/retail.xlsx",skip =1)

## readxl works best with a newer version of the tibble package.
## You currently have tibble v1.4.2.
## Falling back to column name repair from tibble <= v1.4.2.
## Message displays once per session.

myts = ts(retaildata[,"A3349873A"], start = c(1984,4), frequency = 12)
autoplot(myts)

The series has an overall increasing increasing trend and also has a seasonal pattern.

myts %>% seas(x11 = "") %>% autoplot()

The variance in the seasonality is decreasing with time. Also, some outlier can be seen in the remainder as well like in year 1992, 1997, 2003. Other then that the X11 decomposition works quite well for the data.

3. Monthly Canadian gas production

autoplot(cangas) +
  xlab("Year") +
  ggtitle("Monthly Canadian gas production")

ggseasonplot(cangas)

ggsubseriesplot(cangas)

Seasonal pattern with trend can be seen with a constant trend between 1973 to 1987. The gas production amount increased in winter and decreased in summer. Maybe the cold weather in winter increased the demand for the gas that the production amount increased. But the production amount also increased in summer as time went on. Maybe it happended because the demand for electricity to run air conditioners increased over time.

cangas %>% 
  stl(t.window = 13, s.window = 13, robust = TRUE) -> fit
autoplot(fit)

We are selecting s.window manually because the seasonal pattern is changing so much in variation. A smaller value of s.window will it to change more frequently. We can see that the remainder has high variance where the seasonality had a high variance. This shows that the additive decomposition might not be a good choice which stl function uses.

  autoplot(seasadj(fit)) +
  xlab("Year") + ylab("") +
  ggtitle("Canadian gas Production") 

cangas %>% 
  seas(x11="") -> fit
autoplot(fit)

The remainder is smaller than stl function. As we are using the multiplicative method, the seasonality is decresing first then increasing and then decreasing again. we can notice when the variance of the time series is changing little, the additive model is giving us low error but when the seasonality changing variation, then the multiplicative model performing better.

  autoplot(seasadj(fit)) +
  xlab("Year") + ylab("") +
  ggtitle("Canadian gas Production") 

4.Australian quarterly clay brick production.

autoplot(bricksq) +
  xlab("Year") +
  ggtitle("Australian quarterly clay brick production")

The data shows an increasing trend and then a constant trend. Seasonality seems to have constant variation except there are few outlier. We are going to try both the examples.

bricksq %>% 
  stl(t.window = 5, s.window = "periodic", robust = TRUE) -> fit.constant
autoplot(fit.constant)

bricksq %>% 
  stl(t.window = 5, s.window = 5, robust = TRUE) -> fit.change
  autoplot(fit.change)

We can see that when we are using changing seasonality, it also capturing the effect of outlier in 1975.

autoplot(seasadj(fit.constant), series = "Constant Seas") +
  autolayer(seasadj(fit.change), series = "Change Seas") +
  xlab("Year") +
  ggtitle("Australian quarterly clay brick production") +
  scale_color_manual(values = c("Red","Blue"),
                     breaks = c("Constant Seas","Change Seas"))

As can be seen that the constant seasonality results in more peaks and troughs but it also do not pick the effect of intervention variables which should be explained seperately.

fit.change %>% seasadj() %>% naive() %>%
  autoplot()

bricksq %>% stlf(method = "naive") %>% autoplot() +
  xlab("Years")

bricksq %>% stlf(method = "naive") -> fit
checkresiduals(fit)

The residuals are increasing in size with time. Also, they are correlated with each other. Maybe if a method that use bigger window in estimating trend reduces that.

5. Sales of printing and writing paper

autoplot(writing)

The data has an increasing trend with seasonality. Seasonality seems to be constant over time. A drift method work better for estimatinf the trend-cycle component.

writing %>% stlf(method = "rwdrift") %>% autoplot()

6. Sales for a souvenir shop

autoplot(fancy) +
  xlab("year") +
  ggtitle("Sales for a souvenir shop")

The seasonlity seems to increase in size over time. We need to use Box-Cox transformation to control variacne as well.

fancy %>% stlf(method = "rwdrift", lambda = BoxCox.lambda(fancy)) %>% autoplot(PI = FALSE)